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Simulation of a modified Hubbard model with a chemical potential using 
a meron-cluster algorithm 

J.C. Osborn a * 

a Physics Department, University of Utah, Salt Lake City, UT 84112, USA 

We show how a variant of the Hubbard model can be simulated using a meron-cluster algorithm. This provides 
a major improvement in our ability to determine the behavior of these types of models. We also present some 
results that clearly demonstrate the existence of a superconducting state in this model. 



Numerical simulations continue to be the most 
successful way to determine the low energy be- 
havior of many strongly interacting systems from 
first principles. Tremendous advances in simu- 
lation methods have been made over the years. 
However there are still many problems one can 
encounter when performing a simulation. One 
such problem is the effort required to generate 
independent configurations. Many important ob- 
servables are sensitive to the large scale struc- 
ture of a configuration and therefore require re- 
quire an updating scheme that can make signifi- 
cant changes to a large fraction of variables. An 
important example is observables that are sensi- 
tive to topology. Topology cannot be changed by 
small, smooth movements and requires updates 
that in some sense mimic the large scale fluctua- 
tions of the system. 

The autocorrelation time quantifies the num- 
ber of updating steps needed to produce an in- 
dependent configuration. Typically it scales as 
t oc £ z where £ is a relevant correlation length. As 
one approaches a phase transition the correlation 
length diverges so it is important to minimize the 
dynamical exponent z as much as possible. Stan- 
dard update algorithms that change only a few 
variables at a time typically have z around two. 
This can be lowered to z w 1 or lower using im- 
proved updating schemes such as overrelaxation. 
Even with this improvement it is still difficult to 
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make large scale changes. Cluster algorithms are 
designed to change a large fraction of the config- 
uration variables together which greatly reduces 
the autocorrelations and can approach z = 0. 

The idea of updating clusters of variables in a 
Monte Carlo simulation originated from the work 
of Swendsen and Wang for the Potts model |Q. 
They used the cluster representation of the Potts 
model presented earlier by Fortuin and Kasteleyn 
pj to efficiently update large collections of spins. 
This idea was later extended to O(N) spin models 
by Wolff §. 

The reason the algorithm works so well is that 
the cluster variables are a natural representation 
of the long range modes of the system. Thus clus- 
ter algorithms not only provide a means for effi- 
cient simulation but also provide insight to the 
relevant long distance physics. The algorithm for 
classical spin models allows clusters to grow by 
joining neighboring sites based on the alignment 
of the spins. This can generate clusters of arbi- 
trary shapes which are well suited to update these 
models. 

Most of the cluster methods for quantum sys- 
tems are related to the loop algorithm [Q. Here 
the clusters form loops of sites. This structure 
can naturally make large changes to the parti- 
cle worldlines in the path integral. There have 
been many improvements made in the updating 
schemes for loop-clusters || and it continues to 
be an area of active research. 

Despite the continuing advance of loop algo- 
rithms, they arc still limited in the types of mod- 
els that they can simulate. Many models with 
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frustration or fermions develop a sign problem 
when written in terms of loop variables. The sign 
problem arises from large cancellations of config- 
urations with positive and negative weights. In 
frustrated systems of bosons the signs arise due 
to the competing interactions while for fermions 
the signs come directly from the permutation of 
particle worldlincs. In order to extend loop al- 
gorithms to these systems, the clusters must also 
provide a means for reducing these cancellations. 

The meron-cluster algorithm is a method for 
solving the sign problem in some models. It was 
first introduced to simulate the 2-d 0(3) model 
with a 6- vacuum term || . Since then it has been 
extended to many other models including certain 
systems of fermions |7],^| . The fermionic version 
uses the same cluster variables as in a loop al- 
gorithm. The difference is that these clusters 
also encode important information on the fermion 
sign. For certain models we are able to make an 
exact cancellation of all negative configurations 
with positive ones leaving only positive contribu- 
tions to the partition function. 

In this paper we present an overview of some 
recent work on applying the meron-cluster algo- 
rithm to a variant of the Hubbard model. The 
next section contains a brief introduction the 
standard Hubbard model and why it is impor- 
tant to many areas of physics. In section |^ we 
discuss the key concepts of the meron-cluster al- 
gorithm for spinless fermions. A complete expla- 
nation of the method is available in the literature 
. Section || discusses the extension of the algo- 
rithm to fermions with spin. Here we construct 
a variant of the Hubbard model with the same 
symmetries that can be simulated efficiently. In 
section [| we discuss the addition of a chemical 
potential while maintaining the ergodicity of the 
algorithm. Lastly we show numerical results con- 
firming a superconducting phase in our variant of 
the Hubbard model. 

1. HUBBARD MODEL 

The Hubbard model was introduced as a sim- 
ple model for interacting electrons j9j. De- 
spite its simplicity it is believed to exhibit a 
wealth of interesting phenomena of strongly cor- 



related fermions. It has since been used to de- 
scribe a variety of systems including BCS su- 
perconductors, high temperature superconduc- 
tors and superconductor-insulator and metal- 
insulator transitions. 

The Hamiltonian for this model is 

H hub = - 1 E c U s> + c L <v 

<ij> (T=T,i 

- (n;, T +n it i) (1) 

where < ij > stands for all pairs of neighboring 
sites. The creation (cj a ) and annihilation {c ia ) 
operators for a fermion of spin a obey the usual 
anticommutation relations and rii tC! = cj a c i is 
the number operator. We have written the model 
such that it remains half-filled (an average of one 
particle per site) at zero chemical potential (/z). 

The properties of the model depend on the co- 
efficient of the on-site interaction. In the repul- 
sive model (U > 0), the interaction can be viewed 
as the leading order coulomb repulsion. Here the 
model is believed to exhibit properties similar to 
the high temperature superconductors. In partic- 
ular there is the possibility of finding a transition 
from antiferromagnetism to d-wave superconduc- 
tivity. This however has not been confirmed nu- 
merically due to a severe sign problem ]l0| . 

The attractive model (U < 0) serves as a very 
useful effective model for interacting fermions. 
Here the interaction can originate from the 
phonon mediated attraction important to BCS 
superconductivity. It is a very general model and 
can be useful in studying any system of fermions 
with an attraction including neutron matter. The 
attractive model does not suffer from the sign 
problem in traditional simulations as the repul- 
sive model does. It therefore continues to be a 
popular subject of numerical simulation. Despite 
the extensive studies that have been performed, 
a precise verification of a phase transition is still 
elusive. This is mainly due to the effort required 
to simulate at very low temperatures pT|. 
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1.1. Symmetries 

The Hubbard model has two important symme- 
tries that we will want to preserve when we con- 
struct our variant below. The first is the SU(2) S 
spin symmetry defined in terms of the usual spin 
operators S a = J2i with 



sf 



4,1 c ^T 

- (n i>T - n itl ) 



(2) 



When /i = there is also an SU(2) C charge sym- 
metry (also referred to as pseudospin) whose gen- 
erators are given by J a — £^ J" with 
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The operators and Jr create and destroy a 
pair of up and down spin fermions on a site. The 
third component J 3 is related to the particle num- 
ber. The chemical potential couples to this oper- 
ator and explicitly breaks the SU(2) C symmetry 
to the U(l) particle number symmetry. In two di- 
mensions the particle number symmetry can un- 
dergo a Kosterlitz-Thouless transition to a super- 
conducting state. 

2. MERON-CLUSTER ALGORITHM 

The meron-cluster algorithm is explained in de- 
tail in the literature . Here we only present the 
key concepts necessary to illustrate its extension 
to fermions with spin and the inclusion of a chem- 
ical potential. 

2.1. Path integral formulation 

For simplicity we will consider only nearest 
neighbor interactions such that the Hamiltonian 
can be written as H = j> The partition 
function can be formulated as a path integral us- 
ing the Trotter-Suzuki decomposition |l2]] . First 
the partition function is written as the product of 
M transfer matrices 



Figure 1. Non-zero elements of the plaquette 
transfer matrix for neighboring sites i and j be- 
tween timeslices t and t + 



with e = P/M. Next the Hamiltonian is split into 
2 x d parts corresponding to each interaction di- 
rection and whether the interaction goes forward 
from an even or an odd site. The transfer matrix 
can then be approximated as 

exp( —e H ) rj JJexp(-e ^ h i,i+p,) 



JJexp( -e ^2 



i+A 



(5) 



i odd 



Z = Trexp( -p H ) = Tr [ exp( — e H ) ] 



M 



(4) 



Now we introduce an occupation number state 
between each of the transfer matrix operators. 
This gives a total of N t — 2dM "timeslices" each 
containing a set of occupation numbers. The indi- 
vidual nearest neighbor interactions that appear 
between any given neighboring timeslices all com- 
mute and we can thus consider the transfer matrix 
on each neighboring pair of sites separately. 

2.2. Loop-clusters 

We are now interested in the plaquette transfer 
matrix t +i<ninj\ exp(—ehij)\riinj>t for the occu- 
pation number states of neighboring sites i and 
j between timeslices t and t+1. If the interac- 
tion conserves particle number then there are only 
six elements of the transfer matrix that are non- 
zero. These six elements are depicted in Fig. [I]. 
The four circles on each plaquette represent oc- 
cupation numbers with white meaning empty and 
black filled. The two matrix elements where the 
fcrmion hops from one site to the other will have 
a negative sign if it causes a permutation of the 
fcrmion world lines. We will ignore this extra sign 
factor for now and will discuss its effect below. 

We now introduce a set of bond variables that 
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Figure 2. Vertical (upper) and horizontal (lower) 
bond variables and the transfer matrix elements 
they contribute to. The operator corresponding 
to each bond variable is also shown. 



connect certain sites together. The connected 
sites are then constrained to either have the same 
or opposite occupation numbers. A flip of a bond 
then inverts both occupation numbers it touches. 
This way any placement of the bonds can gen- 
erate several occupation number configurations 
through flipping. The probabilities of placing the 
different bond variables is chosen such that the al- 
lowed occupation number configurations are gen- 
erated with the correct weights from the transfer 
matrix. The model can then be viewed as a sys- 
tem of bonds instead of occupation numbers. 

The simplest set of bond variables one can con- 
sider are the two shown in Fig. |[ The vertical 
bonds contribute to the four transfer matrix ele- 
ments shown and the horizontal bonds also have 
four contributions. If we make the approxima- 
tion exp(— ehij) w 1 — ehij then it is easy to see 
what interaction term in the Hamiltonian pro- 
duces these bond variables. The vertical bonds 
simply correspond to the identity term in the 
transfer matrix which is always present. The hor- 
izontal bonds produce the interaction term 



1 

Tlj 

3 2 
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Other possible bonds are examined in |7|jL3|]. Ev- 
ery site is connected to two bonds, one each on the 
forwards and backwards plaquettes in the time 
direction. The bond variables then form closed 



loops of sites which do not change the magnitude 
of the weight of the configuration when flipped. 

2.3. Meron-clusters 

We now need to consider the previously ignored 
fcrmion sign. Whenever we flip a cluster, the new 
state will have the same magnitude of weight as 
the old but might have a different sign. Clusters 
that change the sign of the configuration when 
flipped are called meron-clusters || . The meron- 
clusters provide a mapping between two different 
configurations that cancel and do not contribute 
to the partition function. 

In order to exploit this cancellation we need 
to satisfy two important conditions. First the 
change in sign by flipping any cluster should not 
depend on how the other clusters are flipped. 
This means that we can cancel the contributions 
from meron-clusters independently of the flips of 
all other clusters. Second we want to ensure that 
all negative configurations are canceled with pos- 
itive configurations by flipping a meron-cluster. 
This ensures that the partition function is the 
sum of only positive configurations. 

We now write the partition function as 



Z = J2 W l b ] Sign[b] 



(7) 



lb] 



where W[b] is the weight for a bond configuration 
b and the sum runs over all possible bond config- 
urations. If the conditions presented above for 
canceling meron-clusters hold, then the average 
sign is given by 

= n e (8) 

loops i flips 

with the product extending over all loop-clusters 
for a given set of bonds. The sum over flips is 
simply 



flips 



if mcron — cluster 
2 otherwise 



(9) 



As was already stated, configurations with 
meron-clusters do not contribute to the parti- 
tion function. We can then efficiently simulate 
the partition function by identifying and avoiding 
mcrons. Typically though, the merons will con- 
tribute to observables. We therefore have to allow 
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some to be generated. The observables we will 
consider only get contributions from configura- 
tions with zero, one or two merons. We can then 
reject any update that would create any more. 

Normally the updates favor the generation of 
merons and we would end up with mostly two- 
meron configurations. Therefore we must sup- 
press the merons by reweighting configurations 
with n merons by a factor of p n with p < 1. We 
then tune p so that the simulation spends roughly 
equal time in configurations of each number of 
merons. As the volume gets larger or the tem- 
perature is lowered, the factor p must be made 
smaller and the tuning becomes more difficult. 
However, even for the largest volumes we have 
run on, this is not a significant problem. 

3. FERMIONS WITH SPIN 

The simplest way to include spin into the clus- 
ter algorithm is to have the clusters only connect 
fcrmions of the same spin and then just dupli- 
cate the cluster configuration for one spin to the 
other. If we only use the vertical and horizon- 
tal bonds then the Hamiltonian takes the form 
H = £<«> h>ij,thij,l where fty )(T is the operator 
(§ for the spin a fermions. This can be rewritten 
as @ 

a 

~\~ 2 Si ■ Sj ~\~ 2 J i ' J j 

- 4 v %] v lA v jA v jd (10) 

with riij = rii + rij, Ui = rii,f + n^j and = 
n<i,a — 1/2- The spin and pseudospin operators 
arc defined in (g) and @. 

We will study the model ([!(]) with an on-site 
interaction and chemical potential (|ll|) discussed 
in the next section. The Hamiltonian for this 
model is more complicated than the plain Hub- 
bard model and differs in that it remains strongly 
interacting at U = 0. However it still has the 
same spin and charge symmetries and also ex- 
hibits a similar phase structure including s-wave 
superconductivity as we will show below. The 
main advantage of this model is that it can be 
efficiently simulated with the meron-cluster algo- 
rithm. 



4. EXTERNAL FIELDS 

We now want to add the on-site interaction 

U ( n *.T - [Kid - ~ fiTli . (11) 

The simplest method of including this term is to 
reweight each cluster with 

V Sign(l) = 2exp(-^-7Z75 / )cosh(/9/iWi) 

flips 

± 2eM~US e ) (12) 

where Si is the total number of sites in a clus- 
ter and Wi is the temporal winding. The minus 
sign is taken if the cluster is a meron and plus 
otherwise. For U < we can easily see that this 
weight is always positive and does not introduce 
any sign problem. For U > the weight can be 
negative, however if U > 2/i then the negative 
signs will always cancel in pairs. It is likely that 
in this region the model will be stuck in a charge 
gap keeping it at half-filling. It is possible to go 
to larger \i before the sign problem becomes se- 
vere, but we have not examined this thoroughly. 
We will only consider the attractive model here. 

4.1. Bridge configurations 

The weight ( |l2"l) for a cluster can fluctuate to 
very large values where the simulation can get 
stuck. We therefore need a method to avoid this 
freezing. The method we employ is based on an- 
other application of the meron-cluster algorithm 
— the quantum Heisenberg model in a magnetic 
field. An efficient method for simulating this 
model using meron-clusters was presented in [fl4| . 
Here we will examine the essential features of this 
algorithm and adapt them to our model. 

If we add a magnetic field h in the z direction 
to the standard loop algorithm for the quantum 
Heisenberg model we would then have to reweight 
each cluster with the factor 

exp(/3hWe/2) + exp(-j3hWt/2) (13) 

where Wi is again the temporal winding of the 
cluster. 

In |Q they worked around the problem of 
freezing by including the magnetic field in the x 
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Figure 3. Bond variables and their weights for 
adding a magnetic field in the x direction to the 
quantum Heisenberg model. 



direction. This is equivalent to the original prob- 
lem since the Heisenberg model is symmetric un- 
der rotations. We can insert this interaction on 
every site on the lattice. The transfer matrix for 
this interaction has the weight cosh(/%/2iV t ) if 
the spin does not change between the two times- 
lices and a weight of smh(/3h/2Nt) if it does. 

A set of bond variables that reproduces these 
weights is shown in Fig. [|. The cluster either con- 
tinues as normal, or is broken so that the pieces 
can be flipped independently. However flipping 
individual pieces may introduce minus signs into 
the weight from the antiferromagnetic interaction 
p4| . Pieces of a cluster that start and end on 
sites of opposite spatial parity will change signs 
when flipped. The sign problem that is intro- 
duced is easily solved with a meron-cluster algo- 
rithm. Since the merons are easily identified, one 
can then suppress their creation and avoid a sign 
problem. 

It is easy to see that this method generates 
clusters for the partition function with the same 
weight as (|l3|). Consider a cluster with n sites 
of one spin (e.g. up) and m sites of the opposite 
spin. Now add the weights of all magnetic bond 
configurations that do not produce any merons. 
This means that we can put the broken bonds on 
any combination of even sites or on the odd sites. 
Let w s = exp(—/3h/2N t ) be the weight for a solid 
bond and w& = sinh(/3/i/2iVt) be for a broken 
bond. The total weight is then 

w™(w s + 2w b )" + w^ l (w s + 2w 6 ) m (14) 

where the extra factor of 2 is due to flipping the 
different pieces separately. This is equivalent to 
dl3| ) with Wi = (n — m)/N f . We now look at why 
this method is mode efficient. 



The key feature of this algorithm is the con- 
trolled creation of the merons. If the merons 
were never allowed to form the simulation would 
freeze just like the plain reweighting method. The 
meron configurations then allow the simulation 
to flow much more easily between configurations 
with large weight. The meron configurations do 
not contribute to the partition function and are 
in some sense generated with an "unphysical" 
weight. These types of configurations have been 
referred to as bridges in the literature [ ft5| . This 
naming emphasizes their role in allowing easy 
passage over areas with small weight. Of course 
the merons are important in measuring observ- 
ables and are thus necessary anyway. In this 
example they have an added feature of avoiding 
freezing. 

Just as we summed the weights for the zero- 
mcron sector, we can now see exactly what weight 
the meron configurations are generated with. In 
this algorithm the meron pieces are always pro- 
duced in pairs. The weight of the two-meron con- 
figuration is given by summing over all ways of 
putting the broken bonds such that there are only 
two pieces that have their ends on one odd and 
one even site. This means that we can split the 
cluster into two parts and one part has broken 
bonds only on the odd sites and the other part 
has breaks only on the even sites. Given a split- 
ting of the cluster such that one side has n e even 
sites and the other has m odd sites, the weight 
of that splitting is 

wl + w l - n *- m °(w s + 2w b )"=(w s + 2w 6 ) m » (15) 

where I is the length of the cluster. The full 
weight of the two-meron sector is the sum of ( |l5| ) 
over all splittings while avoiding possible dupli- 
cates. This weight is not as easy to calculate as in 
the zero-meron sector. We can estimate its mag- 
nitude by identifying the leading term for large 
fields. The leading term of ( |l5| ) is 

exp(0h[n e - n a + m - m e ]/2N t ) (16) 

where we have made the substitution I = n e + 
n a + m e + m Q with n the number of odd sites 
on the same side of the splitting as n e and m e 
is the number of even sites on the side with m . 
The terms w n = n e — n Q and w m — m e — m give 
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the temporal extent of each part of the splitting. 
The total winding number of a cluster is given by 
Wg = (w n + w m )/N t and is what appears in the 
zero-meron weight (|l5|). Instead the two-meron 
weight is related to the maximum difference in 
temporal winding of the two parts when the clus- 
ter is split. 

Consider now why the traditional reweighting 
scheme fails. If a cluster with a large positive 
winding is next to a cluster of large negative wind- 
ing the two can never join. This is because the 
cancellation of the positive and negative windings 
would produce a cluster of much smaller weight. 
However if this combined cluster were given the 
weight for the two-meron configuration the weight 
would typically be greater than that of the two 
separate clusters. The clusters would now be al- 
lowed to join. 

Instead of implementing the meron-cluster al- 
gorithm for this model we can just generate con- 
figurations using either the regular (zero-meron) 
weight or the bridge (two-meron) weight. For the 
Hubbard model we need to reweight clusters ac- 
cording to ([[2]) . For simplicity we then choose the 
bridge configuration weight to be 

2exp(-^-i[/^)cosh(/^K - w m ]/N t ) (17) 

with w n and w m defined as above for the splitting 
that maximizes this quantity. This weight has the 
same behavior for large chemical potential as the 
true two-meron configurations would and can be 
calculated efficiently in the update process. 

As with the implementation of the meron- 
cluster algorithm the simulation will want to pro- 
duce mostly bridge configurations. We there- 
fore need to suppress the bridge configurations 
by some constant factor. For large fields choosing 
the right factor becomes difficult, however for the 
results presented below we find good movement 
between the regular and bridge configurations. 

5. NUMERICAL RESULTS 

We simulated the model ( |Io| ) with an chemi- 
cal potential of fi — 1 in two spatial dimensions. 
All results shown are with U = 0. The model is 
still strongly interacting here and has a qualita- 
tively similar behavior as for U < 0. To make 
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Figure 4. Pair susceptibility versus temperature 
for several system sizes. 



sure that the simulations weren't getting stuck 
in any area of phase space we ran at least ten in- 
dependent simulations for each set of parameters. 
Each run used 10 4 configurations for the measure- 
ments. The error bars shown are estimated from 
the variance over the independent runs. 

In two dimensions long range order is forbidden 
according to the Mermin- Wagner theorem. How- 
ever long range correlations can emerge through 
a Kosterlitz-Thouless phenomenon . This can 
be observed by measuring the correlations of the 
pair annihilation operator A = J^. Ci^Ci^ with 
the pair creation operator A^ . The scaling behav- 
ior of the pair correlations is most easily observed 
through the pair susceptibility which we define as 



-(0-m A t e"* H A 



(18) 



In a finite volume for T < T c the pair susceptibil- 
ity should scale with the system size (L) as 



XL oc 



L 7(T) 



(19) 



with 7(T C ) = 7/4 and j(T) approaches 2 as T 
approaches zero. Above T c the susceptibility will 
approach a constant. 

In Fig. U we show the pair susceptibility ver- 
sus temperature for several system sizes. At high 
temperatures the susceptibility does not grow 
with the system size. As the temperature is low- 
ered the susceptibility for large volumes grows 
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o T = 0.2 (y = 1.91) 
o T = 0.526 (y= 1.77) 
□ T = 0.625 (Y= 1.60) 
A T = 0.667 




Figure 5. Pair susceptibility versus system size 
for several temperatures. 
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Figure 6. Scaled winding number versus temper- 
ature for several system sizes. 



rapidly. For low enough temperatures the suscep- 
tibility appears to scale as a power of the volume. 

To see this scaling more clearly, in Fig. |5| we 
show the pair susceptibility versus system size on 
a log-log scale. At the two lowest temperatures in 
the plot (T = 0.2 and T = 0.526) we find excellent 
power law scaling with exponents of 7 = 1.91 and 
7 = 1.77 respectively. The 7 of 1.77 is just above 
the prediction at T c . This suggests that T c is 
just above 0.526. Indeed as we raise the tempera- 
ture to 0.625 we no longer find power law scaling. 
Smaller volumes (L < 32) appear to scale, but 
when we go to larger volumes we see that this is 
clearly not the case. It is therefore crucial to sim- 
ulate on large volumes to accurately determine 
the critical behavior of the system. Conventional 
methods are limited to much smaller lattice sizes 
and therefore have not been able to accurately 
measure any long range pair correlations in the 
attractive Hubbard model. 

Another important observable for studying su- 
perconductivity is the superfluid density. It can 
be defined in terms of the spatial winding as JllJ 

Ps = I < (W x /2f + (W y /2f ) (20) 

where W x (W v ) is the total number of particles 
winding around the boundary in the x (y) direc- 
tion. This quantity is sensitive to the phase co- 
herence of the system and is thus used to indicate 
superconductivity. 



The finite size scaling formula for T < T c is 
given by @ 

^Ps = 2 + ^A coth(v^ log(L/L )) (21) 

with A = at T = T c . For T > T c the superfluid 
density vanishes at large volumes, and jumps to 
a universal value at T c . 

In Fig. H we show the scaled superfluid density 
for different volumes as a function of temperature. 
At large temperatures we can clearly see that it 
goes to zero as the volume is increased. For low 
temperatures we do not see any scaling with the 
volume. At T = 0.526 there is little scaling with 
the volume and a fit to (|Tj) verifies that it ap- 
proaches a constant at infinite volume. This is 
consistent with the conclusion reached from the 
pair susceptibility that this temperature is below 
T c . However at T = 0.555 there is a noticeable 
scaling and the fit predicts that the superfluid 
density should scale to zero. 

To show that we are simulating a fermionic 
model we can look at the double occupancy which 
is twice the average number of sites that have 
both an up and down spin fermion. If the 
fermions are all tightly bound in on-site pairs then 
this quantity would be equal to the total parti- 
cle number. In Fig. we see that the average 
number of particles stays roughly constant just 
above 1.1 as the temperature is changed. The 
double occupancy, however, starts off fairly low 
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versus temperature. 



at high temperatures and rapidly increases as T c 
is approached. Thus we are able to observe the 
formation of the pairs responsible for supercon- 
ductivity and not just simulate the Bose-Einstein 
condensation of preformed boson pairs. 



6. CONCLUSIONS 

We have shown how a meron-cluster algorithm 
can be used to efficiently simulate a variant of 
the Hubbard model on large lattices. This has al- 
lowed us to accurately study the behavior of this 
model near the superconducting phase transition. 
Currently the meron-cluster algorithm is limited 
to specific models. However, just as cluster algo- 
rithms for bosonic systems have flourished over 
the past ten years, hopefully fermionic cluster al- 
gorithms will too. 

We are currently using the model shown here 
to examine the superconductor-insulator transi- 
tion in the presence of disorder. Preliminary re- 
sults suggest the possibility of a localized super- 
conducting state, though further work is needed 
to determine the precise nature of this state. An- 
other future interest is extending the fermion con- 
tent to include isospin or color degrees of freedom. 
This may provide some precise results for new and 
interesting types of models. 
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